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ABSTRACT 

Ever since the discovery of the edge-on circumstehar disk around /3 Pictoris, a standing question 
has been why the gas observed against the star in absorption is not rapidly expelled by the strong 
radiation pressure from the star. A solution to the puzzle has been suggested to be that the neutral 
elements that experience the radiation force also are rapidly ionized, and so are only able to accelerate 
to an average limiting velocity ?;ion- Once ionized, the elements are rapidly braked by C II, which is 
observed to be at least 20 x overabundant in the disk with respect to other species. A prediction from 
this scenario is that different neutral elements should reach different Wion, depending on the ionization 
thresholds and strengths of driving line transitions. In particular, neutral Fe and Na are predicted 
to reach the radial velocities 0.5 and 3.3kms~^, respectively, before being ionized. In this paper we 
study the absorption profiles of Fe and Na from the circumstehar gas disk around (3 Pic, as obtained 
by HARPS at the ESO 3.6 m telescope. We find that the Fe and Na velocity profiles are indeed shifted 
with respect to each other, confirming the model. The absence of an extended blue wing in the profile 
of Na, however, indicates that there must be some additional braking on the neutrals. We explore 
the possibility that the ion gas (dominated by C II) can brake the neutrals, and conclude that about 
2-5 X more C than previously estimated is needed for the predicted line profile to be consistent with 
the observed one. 

Subject headings: circumstehar matter - protoplanetary disks - stars: formation - stars: individual 
{J3 Pictoris) 



1. INTRODUCTION 

Understanding the evolution of gas in circumstehar 
(CS) disks is critical to understand star and planet for- 
mation. Unfortunately, CS gas in general is very diffi- 
cult to observe, so how and when the gas disappears, 
is still poorly known. The main difficulty is that the 
cool gas in disks only thermally excites the lowest ly- 
ing energy levels, which emit radiation in mid- to far- 
infrared wavelengths that are blocked by the terrestrial 
atmosphere. Radiatively excited levels are easily ob- 
served from the ground, but require very high contrast 
sensitivity, as the star typically emits many orders of 
magnitude more radiation at the same wavelengths. In 
the case the CS disk is seen edge on, however, the gas 
can be seen in absorption against the star. This is 
the case for the disk arou nd Pictoris, w here CS ab- 
sorption lines were foun d flSlettebaklll975D long before 
the d isk was discovered (jAuman nlll985l: ISmith fc Terrild 
Tp3). /3Pic is a nearby (19.44±0.05 pc: Ivan LeeuwenI 



2007) A5 main -sequeiice sta r estimated to be 10-2 Myr 
old (Zuckerm an et all [20011 : iMentuch et al.l |2008f). Be- 
cause time scales for dust removal mechanisms are much 
short er than the age of the system (jBackman fc Parescd 
|1993[ ). the dust disk must be replenished. The most pop- 
ular mechanism to replenish the dust disk is through col- 
lisional cascades from planetesimal size objects, hence 
the classification of the (3 Pic disk as a debris disk. 

Even though the fate of CS gas is much less con- 
strained than the dust, also the gas disk around (3 Pic 
is thought to be secondary as opposed to a remnant 
from the star forming nebula, as indicated by the pres- 

* Based on observations made with ESO Telescopes at the La 
Silla Observatory. 



ence of C O (iVid al-Madi ar et al.lll994f) that quickly is dis- 
sociated (Ivan Dishocck fc Blackl 1988f). and bv dv nami- 
cal a rguments ([Fernandez. Brandeker. fc Wull2006l here- 
after |FBW06[ ) . Suggested mechanisms f or replen ishing 
the gas are photo-desor ption of dust (Ch en et al.l [2b07). 
evaporation of comets (Bcust et al. 1989), and collisional 
evaporation of dust grains (Czechowski fc Mann 2007). 

After the discovery of the /3Pic CS dust disk, more 
attention was directed to wards the gas obs^ erved in ab- 
sorption against the star ([Hobbs et al.lll985l) . and it was 
soon realized that the expelling radiation force on many 
elements exceed the gravity — by up to a factor of several 
hundred in the case of Na I. Yet, most of the gas absorp- 
tion is observed to be at rest with respect to the star. 
Evidently, there must be some braking agent at work, 
but the lack of spatial information of the location of the 
gas made it difficult t o constrain braking models — e.g. 
ILagrange et al.l ([1998 ^) suggested that a dense ring of H I 
at sub-AU distance from the star would be able to brake 
the gas through friction, without violating any observa- 
tional constraints known at the time. This changed with 
the detection of spatially resolved Na I emission from the 
disk, observed to be co-located with the dust disk out to 
hundreds of AU (Olofsson et al. 20 01j). For H to b rake 
the Na I over the whole disk, IBrandeker et al.l (|2004[ ) es- 
timated that the required mass woul d have to be h igher 
tha n upper limits set on H I (Frcudli ng et al.|[l995D and 
H2 ([Lecavelier des Etangs et al.i i2001). 

The problem of how the gas arou nd (3 Pic is braked 
was investigated in detail bv IFB W^ . One key observa- 
tion they made was that all neutrals that are subject to 
a strong radiation force, are also quickly ionized. The 
neutrals thus only reach small velocities Vion on average 
(depending on element) before being ionized. Since ions 
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are affected by long-range Coulomb forces, they interact 
strongly and act as a single fluid; ions experiencing radi- 
ation force are efficiently braked by those unaffected by 
radiation. For ions, C II turns out to be the most dom- 
inant braking agent, while Fe II is the strongest driver. 
For gas of cosmic abundance, the effective radiation force 
acting on the g whole, was found to be 4 x stronger 

than gravity. In effect, Fe II would drag the rest of the 
gas along and expel it from the system. In order for the 
ion gas to be stable against the radiation pressure, an 
overabundance of at least 10 x of C would be required, 
and in fact, FUSE observations of C in absorption against 
the star subsequently proved C to be at least 20 x over- 
abun dant, providing the stabilizing inertia to brake the 
disk (jRoberge et al.li2006l hereafter EOBOl). 

A prediction from the scenario that only ions are 
braked, is that neutrals should be observed to reach a ra- 
dial velocity of uion with respect to the star. The purpose 
of this paper is to test that prediction in detail, by mak- 
ing use of archival high spectral resolution data that has 
been obtained of /? Pic for other purposes. Since the abso- 
lute heliocentric velocity of /3 Pic is very hard to measure 
to good accuracy we treat it as a free parameter and fo- 
cus on the observed differential velocity between Fe I and 
Na I. The paper is organized as follows. In §[5] we study 
what the implications of the braking scenario look like 
in detail, in particular the expected absorption line pro- 
files. §[3]reports on how the data we use to compare with 
theory were acquired and handled, especially the some- 
what elaborate procedure to remove telluric lines from 
the Na I region. In §[5] we show the resulting compar- 
isons, and discuss implications. The paper is concluded 
by §1 

2. LINE PROFILES 

If disk braking operates on ions, but not on neutrals, 
we can expect the neutrals to be accelerated for a short 
while before getting ionized. Since the number of ac- 
celerating photons outnumber the ionizing photons by 
at least two orders of magnitude, we can treat accel- 
eration as a continuous process while the time before 
ionization is stochastic according to an exponential dis- 
tribution. Given an ensemble of particles, this will give 
rise to a distribution of velocities at any given moment, 
which will broaden an absorption line and result in a 
characteristic profile. In the following two subsections 
we first study the case that there is no braking acting 
on the neutral particle, and then continue with the case 
where some braking is provided by a surrounding ion gas. 

2.1. Freely accelerated neutrals 

To find out what kind of absorption profile we should 
expect from neutral elements being radially accelerated, 
we make the following assumptions: 

1. The neutral particles are ionized at a rate T(r) cx 
r~^, where r is the distance to the star. The par- 
ticle has no "memory" ; the lifetime A of a neutral 
particle is thus exponentially distributed with the 
expected lifetime F"^. 

2. The particles are in ionization equilibrium, mean- 
ing that the recombination rate equals the ioniza- 
tion rate. 



3. The particles travel a very short distance before 
being ionized; r is thus essentially constant during 
the neutral phase. 

4. Ionized particles are braked efficiently and so are 
at zero radial velocity. Neutral particles thus start 
with no radial velocity. 

The probability density distribution of the lifetime A be- 
fore ionization is thus 



/A(t) = rexp(-ri). 



(1) 



The radial equation of motion for a particle when it be- 
comes neutral is 



1) 



(2) 



where /? = -Frad/-Fgrav is the ratio between radiation force 
Frad and gravity Fgrav, G is the gravitational constant, 
and Mi, is the mass of the central star. The angular mo- 
mentum L = (GAfro)^/^, assuming the particle becomes 
neutral while in a circular Keplerian orbit of radius tq. 
For a particle moving a short radial distance before be- 
coming ionized again (and thus r « ro), the radial accel- 
eration during its neutral lifetime can be approximated 
with 

a = (3—^, (3) 

The velocity reached by a particle before ionization is 
then Vion = oA, and its velocity probability density dis- 
tribution 

.-1 ' ^ 



(4) 



where we have substituted v = at and Uion = a/P is the 
expected radial velocity of a particle before ionization 
(and so the radial dependence r~^ is divided out here). 

Since the lifetimes are exponentially distributed, they 
do not depend on when the estimate is made, the dis- 
tribution of remaining lifetime is always the same. In 
particular, since we have ionization equilibrium, the life- 
time probability distribution is equal to the distribution 
of times since recombination. The normalized column 
density profile of particle density as a function of veloc- 
ity N{v) is therefore equal to Eq. UJ 



N{v) 



, exp(- 



-)• 



(5) 



2.2. Neutrals braked by ion gas 

If the neutrals moving radially outwards are interact- 
ing with a gas at rest, the neutrals will be braked. The 
neutral-neutral interaction cross-section is tiny, however, 
and requires a very high gas density to be effective. Much 
more efficient is neutral-ion interaction, due to the po- 
larizability of neutral elements. Since the gas around 
/3 Pic is strongly ionized, we now look closer on how the 
absorption profile would be affected from the braking. 

The braking of neutrals in an ion gas is analogous to 
ions being braked in a neutral g problem studied by 
iBeust et al.l ()1989[ ). The equation of motion for a neutral 
particle of mass m accelerated by radiation and braked 
by an ion gas is 



dv ^GM^^m , 

m— = p kv, 

at 



(6) 



Gas braking in the /3 Pic disk 
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where 



/Mi 



(7) 



Ui is the number density of ions of species i, rrii is 
the ion mass, = mim/ [rrii + m) is the reduced mass, 
a is the polarizabihty of the neutral particle, e is the 
electron charge, and eo is the electrical permettivity. The 
asymptotic stable solution is 



If Wdrift^drag < T, whcrC ^drag = m/k, 

dependent solution is 



(8) 

then the time- 



V{t) « Wdrift [1 - exp(-i/idrag)] 



(9) 



where we have assumed that the particle starts without 
initial velocity. The probability distribution for the life- 
time A of a neutral particle is still given by Eq. [TJ so the 
resulting velocity distribution f (A) is described by the 
probability density distribution 



N{v) - /a {t{v)) 



dt{v) 



Av 



1 - 



(7-1) 



(10) 

where 7 = Wdrift/^ion — Ttdrag- For 7 > 1, Eg.fTUlreduces 
to Eq. [S] The expected velocity is 



which turns to 



E[v]^ 



E[v] 



IVic 



7+1 



(11) 



(12) 



for 7 ^ 1 . From Eq. [TO] follows that there generally is a 
radial dependence of the velocity profile (as it depends on 
local circumstances, such as ion gas density), in contrast 
to the freely accelerating profile of Eq. [5] 

2.3. Self- shielding: numerical model 

The photons that contribute to the acceleration of the 
gas are generally scattered in a direction different from 
the line of sight. This extinction implies that gas lo- 
cated further out will see fewer photons, resulting in a 
lower effective radiation force coefficient /3; this is called 
self- shielding. The coefficient /3 for a particle will thus 
depend on the column density of particles between its 
position and the star, of particles that share the same 
velocity. The /? coefficient then in turn determines the 
acceleration, and what velocity the particle will reach 
before getting ionized. Because /3 is changing with ve- 
locity, and depends on the radial density profile of par- 
ticles, we model this effect of self-shielding numerically. 
The numerical model is simplified by the one-dimensional 
nature of the problem, and that gas at outer radii are 
affected by gas inside, but not the reverse. For mod- 
eling the absorption profile from /? Pic, we assume the 
density, ionization, and temperature distribu tion com- 
puted bylZagorovskv. Brandeker. fc Wul (|2010l hereafter 



putea by 
IZBWld l. 



The procedure adopted as follows: 



1. Divide the radial distance into A^shcU shells of thick- 
ness Ar. Start with the innermost shell and com- 
pute the velocity distribution of particles in each 
shell going outwards. 



2. Use A'^part number of particles in a shell (where 
each particle represents a number of real parti- 
cles). Initialize all particles with a random ther- 
mal velocity v according to the probability dis- 
tribution Vt{v) = (27rcr|)-i/2exp [-0.5(i;/crT)2], 

1/2 

where ctt = [AjBli/m] , is the temperature of 
the shell i, and m is the mass of the particle. The 
radius of shell i is r^. 

3. Use the velocity dependent optical depth t{v) to- 
wards the star to compute /3 as a function of v. 
There is no extinction inside the first shell. 

4. Iterate Mtcr times. For each iteration, evolve the 
velocity of each particle according to the equation 
of motion of Eq. [S] using a time step of At oc 
I'ion with the addition that /3 is a function of 
V. For each iteration and particle, reset the parti- 
cle velocity to the thermal velocity with probability 
FionAt. This simulates that the particle gets ion- 
ized. Since we assume ionization equilibrium, there 
is a recombination for each ionization, and the re- 
combined particle is set to the thermal velocity of 
the gas (the assumption is that ions are efficiently 
braked to on average ^zero radial velocity). Af- 
ter Alitor iterations the solution has converged to a 
stable statistical distributions of radial velocities of 
the particles. 

5. Make a histogram of the velocities of all particles 
inside the present shell, with A'voi velocity bins and 
bin size Aw. Normalize the distribution such that 
it corresponds to the total column density of the 
sheU, ^N,{vj) = ngas(r,)Ar. 

6. Compute the extinction towards the star for the 
next shell by adding the velocity dependent column 
density for each shell, Ni{vj) = J2k<i '^^kivj), 
and find Ti{vj) = Ni{vj)aoc/{XoAv), where Aq is 
the wavelength of the transition and 



Aij Aq ffj 

Sttc gi 



(13) 



is the cross-section times wavelength per particle 
for the transition, where Aij is the Einstein coef- 
ficient and gj and gt the statistical weights of the 
levels. 

7. Repeat from ([2]) for all shells. The observed profile 
can then be estimated from Tjv^j^^j|(t;). 

We implemented this model in matlab for Na I around 
/3 Pic, and verified it by testing against conditions where 
the analytical profiles of §S I2.1I and 12.21 should be valid, 
with excellent agreement. In detail, we used A'shcii = 186 
with n between 15 and 200 AU (and thus Ar = 1 AU), 
A^part = 10'' per shell, A^iter = 10^ iterations, and a time 
step of At = [2000rion(r)]~"' (meaning we iterate 5 times 
longer than the ionization time scale). For the velocity 
grid, we used Av = 100 ms~^ from —2 to 20kms~^, 
to be certain to cover most of the particle velocities. 
For the Na I D2 line ao = 1-987 x 10"" cm^ atom'^ 
and for the Di line the cross section is 1/2 as much. 
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TABLE 1 
Lines studied 



Element 


Design. 


Aair[A] 


A„,[s-i] 


9u 


-9l 


Ref.': 


Ca II 


K 


3933.6614 


1.47 X 10* 


4 


- 2 


1 


Ca II 


H 


3968.4673 


1.4 X 10** 


2 


- 2 


1 


Fe I 




3824.4436 


2.83 X 10^ 


7 


- 9 


1 


Fe I 




3859.9114 


9.70 X 10^ 


9 


- 9 


1 


Nal 


Da 


5889.950954 


6.16 X 10^ 


4 


- 2 


2 


Na I 


Di 


5895.924237 


6.14 X 10^ 


2 


- 2 


2 



Wav elength references are (1) IMortonI 1(2003 ). (2) IWolf et aLl 
((2003) 

Since these two resonance lines are completely domi- 
nating the radiation force on Na I, we computed the 
optical depth T{vj) for the D2 line and the effective 
Pivj) = /3Nai(exp[-T(uj)]2/3 + exp[-r(wj)/2]/3), where 
/^Nai is the radiation force coefficient on Na I without 
extinction. 

3. HARPS DATA 

To test how well the ion braking scenario compares 
to observations, we retrieved 489 spectra of /3Pic from 
the ESO archive. The spectra were obtained by HARPS 
at the ESO 3.6 m telescope at La Silla, Chile, dur- 
ing the period of 2003-10-28 to 2008-03-21. HARPS 
is a fiber- fed high-resolution [R ^ 10^) echelle spec- 
trograph (A = 3800 — 6900 A), dedicated to the search 
for exoplanets. The long-term radial velo city accuracy 
of HA RPS is estimated to be < Ims""'^ (jMavor et al.l 
I2003f) . The HARPS archive at ESO contains data pro- 
cessed by the HARPS science-grade reduction pipeline, 
producing high-quality extracted wavelength-calibrated 
spectra. Our main interest are regions around the Ca II 
H & K doublet at 3934 A and 3968 A, the Fe I 3860 A 
line, and the region around the Na I D doublet at 5890 A 
and 5896 A (see Table[l]for line details). Typical signal 
to noise ratios in single spectra range from 50-100 for 
the 3900 A region, to 100-250 around 5900 A, depending 
on exposure time and seeing. 

3.1. Cleaning telluric lines 

With the high-quality reduced data produced by the 
HARPS pipeline, not much post processing needs to be 
done, in general. The exception is the Na I absorption 
region at 5890 A, where abundant time- varying absorp- 
tion features from atmospheric H2O interfere with the 
CS absorption lines of interest. Common practice is to 
remove the telluric features by observing a nearby "tel- 
luric standard" right before or after the sc ience obser- 
vatio n, and divide out the features (e.g., iVacca et al.l 
120031) . Another is to fit a theoretical atmosphere spec- 
trum ()Seifahrt et al.|[2010[ ). Since we did not have access 
to a standard star and it is hard to make a good model 
fit without leaving systematic residuals, we opted for a 
third method, making use of the strengths of this partic- 
ular data set: 

1. The supreme stability of HARPS (including the 
line-spread function). 

2. The large number of observations (489). 



3. The fact that all telluric absorption lines in the 
region are optically thin. 

4. The variation in position and strength of the tel- 
luric features, due to the variable barycentric cor- 
rection from the orbital and rotation speed of the 
Earth and the variable airmass and water content 
of the atmosphere. 

By assuming that the stellar spectrum is static and that 
only the strength and relative shift of the telluric lines 
are varying, we can use the large number of spectra to 
constrain the problem. To constrain the telluric absorp- 
tion lines even further, we downloaded multi-epoch spec- 
tra for three additional AO stars: HD 71155 (40 spectra), 
HD 177724 (12), and HD 188228 (24). The advantage of 
AO stars is that they have very weak Na I features that 
are rotationally broadened to several 100kms~^, making 
them suitable as spectrally smooth background lamps to 
study telluric absorption lines. The problem now be- 
comes to find a (normalized) telluric extinction spectrum 
k(A), a model spectrum F^^^{X) (for each star fc), and a 
water column parameter ak^n (for each observation k, n) 
such that 

Xmod = X! X! / fc.n/,N l-^obs (^) - 

^mod([l + ^fe,n]A) exp(-afc,„K[A])|dA (14) 

is minimized, where Nk is the number of observations of 
star k, is the error of the measurement of spec- 

trum k, n at wavelength A, and Zk,n is the (known) red- 
shift associated with observation k, n. We implemented 
a non-linear conjugate gradient minimization method in 
C-I--I-, using the Polak-Ribiere method with res tart to 
guarantee convergence (see, e.g.. lShewchukl[T99^ . Since 
minimization in general is degenerate with respect to 
a continuous spectrum (i.e. many different spectra can 
produce the same minimum x^), we added a second- 
order "smoothness" regularization that penalizes rapid 
changes of the slope: 

X?eg({/™}) = E (l^I^^^±I^ -fS, (15) 

where fm is the value of pixel m of the discrete model 
spectrum (or model extinction) of M pixels, and fo — 
/m-1 are boundary values (arbitrarily fixed to 1 for stel- 
lar spectra and for extinction). In practice we thus 
minimize 

= Xmod + aXrcgill^m}) + bJ2Xrcgi{Pmod,m}), (16) 

fc 

where a and b are regularization constants that decide 
the relative importance of the regularization. These con- 
stants have to be manually tuned to produce good re- 
sults; too small and the minimization process will pro- 
duce numerically noisy spectra, too high and the spec- 
tra will be too smooth, devoid of detail. Fortunately, 
there is a large range between the too small and the 
too high, that produce good spectra that only weakly 
depend on the regularization constants. We found that 
a = b = (10 — 1000) X -/Vgpec, where A^spec is the total 
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Fig. 1. — The upper panel shows a single spectrum of /9Pic as 
obtained from the HARPS pipeline around the Na I D lines. This 
region is plagued by numerous telluric absorption lines from H2O. 
The lower panel shows the resulting spectrum from a combination 
of 92 FEB-fr ee sp ectra, where the tellurics have been cleaned (as 
described in ^ 13. 11 1. 

number of spectra used, produced good results; in the 
end, we used a = b = 250 x iVspec- 

After an optimal fit to the telluric lines were found 
(simultaneously with the average stellar spectra), we di- 
vided ( "cleaned" ) all individual observations by the fitted 
telluric extinction. This allows us to inspect individual 
clean spectra and accommodate potential changes be- 
tween epochs. The improvement from the cleaning pro- 
cess is illustrated in Fig. [1] 

3.2. Deriving line profiles 

In addition to the stable gas absorption lines 
seen against /3Pic, some lines are known to show 
strongly velocity-shift ed transient absorption features 
([Lagrange et al.( [l987D . These absorption features are 
thought to be originating fror n "falling evapora ting bod- 
ies" (FEBs) close to the star (Beust et al."1989'), and are 
mostly observed in the UV (Dclcuil ct al. 1993) and in 
the Ca II H and K lines. Their relevance to the present 
work is that we want to avoid spectra with strong FEB 
features, as they may bias the line profiles we are study- 
ing. Fortunately, neither Na I nor Fe I are known to 
show FEB features, but that is likely a sensitivity issue. 
Since Ca , Na, and Fe around f3 Pic are strongly ionized 
(jZBWlOl ). the column density of Ca II is much higher 
than either Na I or Fe I, making the Ca II H and K 
lines strongly saturated and FEB features in those lines 
stronger, despite having oscillator strengths similar to 
the Na I and Fe I resonant transitions. To reduce the pos- 
sibility of FEBs affecting the Na I and Fe I line profiles, 
we inspected all 489 spectra for FEB features in the sensi- 
tive Ca II lines, and selected 95 without any sign of FEB 
activitjQ. With no FEBs visible in the Ca II lines, we 

^ In order to provoke a FEB signal in the Na I lines, we also se- 
lected a sample of 262 spectra that show strong red-shifted FEB ac- 
tivity in Ca II. Adding the frames, the resulting Na I indeed shows 
a red-shifted FEB feature with an equivalent width O.S^q^ mA, 
i.e. at the ~10% level of the main stable line. 
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Fig. 2. — The observed normalized line profiles of Na I Di and 
Fe I A3860 with over plotted models. The crosses show the observed 
Fe I line profile, the plus signs the observed Na I line, and the 
solid lines are from the nominal model (described in §|4]l. As seen, 
the blue wing of the Na I line is greatly overestimated by the 
model compared to the observed profile, while the Fe I line looks 
fine. The dashed line shows the best-fit model for Na I, where the 
carbon abundance has been increased by a factor of 2 compared 
to the nominal case (the Fe I profile does not change; see JSj. The 
expected errors on the measured profi les a re 1 % of the peak for 
Na I and 3 % of the peak for Fe I (see S l3^ . 

expect their contribution to the Na I and Fe I lines to be 
negligible. Subregions of the 95 spectra around the Na I 
5890 A, Fe I 3824 A, and Fe I 3860 A lines were then ex- 
tracted and interpolated to a velocity scale of 0.1 kms~^ 
resolution (using the HARPS wavelength calibration and 
line center wavelengths from Table [T]), stacked (weighted 
by their signal to noise), and finally normalized. The re- 
sulting profiles for Fe I 3860 A and Na I 5890 A are shown 
in Fig.[2l As seen in the figure, there is a significant offset 
between the observed radial velocities of Na I and Fe I, 
as expected from the ion-braking scenario. 

Given the large number of high signal-to-noise spec- 
tra being combined, errors in the derived profiles are ex- 
pected to be dominated by systematic errors, such as 
flat-fielding, flux calibration, and telluric line modeling 
errors. To check the reliability of the derived profiles, 
we compared two lines from the same element with each 
other; many of the expected systematic effects should 
change with wavelength. One notable exception is con- 
tamination in the line profile due to FEB activity; since 
only spectra without FEBs visible in the Ca II are se- 
lected, we estimate the impact on the Fe I/Na I regions 
to be < 1%. By comparing the Fe I 3824 A with the 
3860 A line, we find them to be consistent to within 3 % 
(1 standard deviation) of the peak. Making the same 
comparison between the Na I 5890 A and 5896 A lines, 
we find them to be within 1 % of each other. The rel- 
ative radial velocity error is negligible in comparison to 
the line profile; it is dominated by the 10 ms^^ accuracy 
by which the Fe I lines are known (HARPS instrumental 
errors are < lms~^). This means that the velocity dif- 
ference between the Fe I and Na I lines, which from Fig. [5] 
is seen to be on the order of 1 kms~^, is lOOx larger than 
the estimated velocity uncertainty, and therefore highly 
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Fig. 3. — The central part of the Na I Di line. The curves 
over plotted the data are the rotational+lS+CS profile, and the 
rotat iona l+IS profile (used to fit the continuum for the CS line). 
See §[33] 

significant. 

3.3. Removing interstellar absorption 

The narrow CS absorption lines in Fig. [1] are seen to 
have "shoulders" , which a re stable in time and due t o in- 
terstellar (IS) extinction (jVidal-Madiar et al.lHIsl ). To 
remove them from the profile, we assumed they were 
Gaussian shaped, and simultaneously fit a rotational pro- 
file with two Gaussians per Na line, one for the the CS 
and IS contribution, respectively (Fig. [3]). The rotational 
profile was best fit by v sin i = 126 ± 1 km s""'^ and a limb- 
darkening coefficient e = 0.57 ± 0.06. The ratio between 
the two IS D2 and Di line strengths is 2.0, indicating op- 
tically thin absorption. The IS lines are centered at a he- 
liocentric radial velocity of 19.46 ± 0.27km s~^ (i.e. close 
to the stellar velocity at 20.5 km s~^) and are unusually 
broad with a broadening parameter b = 21.4±1.7kms~^, 
related to the standard devia tion a of the norma l distri- 
bution through b = iVidal-Madiar et all Hl986f ) 
found the IS line to be present also in the nearby star 
a Pic and were able to remove it from their /3 Pic spec- 
trum (along with telluric lines) by dividing the (3 Pic 
spectrum with that from a Pic. Since we do not have 
access to an a Pic spectrum of a quality similar to the 
/3 Pic spectra, we remove the IS line by division with the 
fitted Gaussian. 

3.4. Estimating the line-spread function 

To compare the theoretical line profile models with the 
data, we need to convolve the models with the instrumen- 
tal profile of the spectrograph, the line-spread function 
(LSF). One way to estimate the LSF is to measure the 
shapes of spectral lines from a Th-Ar lamp, since they 
are expected to be unresolved. Because HARPS records 
Th-Ar spectra in a parallel fiber with (almost) all sci- 
ence observations, we have a huge data set of spectral 
lines to choose from. Since we are interested in spectral 
regions around 3900 A and 5900 A, we decided to deter- 
mine the LSF for those two regions separately by only 
selecting suitable lines (unblended lines of sufficient flux 




Fig. 4. — Using a model where the ionization velocity t)ion is 
reduced as much as possible while still being barely consistent with 
observations, we here plot the square sum of the residuals between 
a fit of a particular 7 as a function of 7/70, with 70 being the 
nominal value. The corresponding column densities are given on 
the upper horizontal axis. 

that do not come close to saturate) from the respective 
orders. With over 500 spectra retrieved from the archive, 
we ended up with several thousand lines for each region. 
The lines were then linearly interpolated and shifted to 
a grid of 0.1 kms~^ resolution, using the wavelength cal- 
ibration provided by the HARPS pipeline. Finally, the 
lines were added, weighted by the signal to noise, fit by a 
third-degree polynomial (for the background) -I- a double 
Gaussian of the functional form 



LSF(u) = fc,exp 



mi 



i=l,2 



2af 



(17) 



where v is the velocity shift from the center of the 
line, and the other constants are fitted. Comparing 
the LSFs determined for the 3900 A and 5900 A regions, 
we found no significant difference (<1% of the peak). 
The fitted parameters ki = 1.193, mi = 5.64ms~^, 
(71 = 1.135kms"\ k2 = -0.194, TO2 = -9.01ms-\ and 
(72 — 0.680 km s^^ approximate the LSF very well, with 
residuals again below 1% of the peak. The full- width half 
maximum (FWHM) of the LSF is 2.67 kms^^, meaning 
a spectral resolution of i? ~ 112 000. 

4. ANALYTIC METHODS 

With a stellar mass of Ah = 1.75 M0 (ICrifo et al.l 
[T99l . [FBW06I compute /^Nai = 360 ±20, /3fcI = 27±2, 
rNai(lOOAU) = 1.1 X 10-^s-i and rFci(lOOAU) = 
5.8 X 10-8 s-i, resuhing in vfj^^ = 3.3 km s^^ and = 
0.5kms~^. As our nominal gas disk model, wc use the 
/3 Pic model with parameters described in ZBWIO. As- 
suming free acceleration, the associated line profiles of 
Eq.[5] are adjusted for the radial velocity of /? Pic by fit- 
ting the Fe I profile (as seen in Fig. [2j the fitted radial 
velocity of the system is 20.5 ± 0.2kms~^). The the- 
oretical profile looks fine for Fe I and the red wing of 
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Na I, but the modeled Na I blue wing is clearly more ex- 
tended than the observed. Apparently the Na atoms of 
the /3 Pic disk do not reach as high speeds as we expect 
before getting ionized. The discrepancy between model 
and observations could be due to the following three rea- 
sons: 

1. Less acceleration, e.g., an overestimated /3 value 
for Na I. Since the j3 Pic spectrum is very well con- 
strained around the Na I D lines, however, and the 
corresponding radiative transition coefficients are 
known to very high accuracy for these resonance 
lines, it is unlikely that the /3 value is more off than 
the estimated 8%. 

2. Higher ionization rate of Na I. This will give less 
time for the atom to accelerate before getting ion- 
ized. The ionization rate is directly proportional to 
the ionizing UV radiation from the star. The UV 
part of the /3 Pic spectrum is not well constrained 
outside the spectral windows observed by FUSE 
(905-1187 A) and HST/STIS (1459-2888 A), but 
most of the ionizing photons are expected to arise 
from close to the ionization threshold at 2413 A, 
where the spectrum is known to within 10% from 
STIS (Aki Ro berge, pri v. comm. 2010). The ioniza- 
tion model of IZBWIOI reproduce observed ioniza- 
tion ratios in the /? Pic gas disk reasonably well, so 
a much different ionization rate than the assumed 
would be surprising. 

3. More efficient braking. With stronger braking, the 
resulting drag velocity is smaller. The braking is 
the strongest for the higher speed particles, possi- 
bly explaining the observed absence of a strong blue 
wing. Since the braking is dominated by, and di- 
rectly proportional to, the number density of C II, 
an underestimate of C II would result in a propor- 
tional underestimate of the braking force. Strongly 
saturat ed abso rption lines of C II observed by 
FUSE (jROBOQl have determined the C abundance 
in the [3 Vic disk to be at least 20 x overabundant 
with respect to other elements in the disk at cosmic 
abundance. Since the stable CS lines are saturated 
and blended with a time- variable broad component 
due to FEBs, small variations in the intrinsic line 
shape can mimic large abundance changes. It thus 
seems motivated to explore the possibility that the 
carbon is even more overabundant than previously 
thought, and what consequences that has for the 
braking and the line profile of Na I. 

A combination of the above is also possible. In partic- 
ular, should braking be more efficient than anticipated, 
the lower radial velocities will contribute to stronger self- 
shielding and a decrease in the /3- value, resulting in less 
acceleration. 

To investigate the possible scenarios producing the ob- 
served line profiles, we made an extensive parameter 
search for Vion and 7, trying to fit the profile from Eg. 1101 
to the observations (leaving the system radial velocity as 
a free parameter). Because self-shielding, as explored in 
the numerical model, in most of the parameter space is 
very small (with an optical depth less than 0.1 resulting 
in only a minor correction to the analytical profile), we 



only used the numerical model to investigate a subset 
of models where self-shielding was expected to be most 
important. 

5. RESULTS AND DISCUSSION 

Since the expected radial velocity increases with both 
Wion and 7, and since the observed line profiles are not 
well resolved, the problem of finding the optimal Wion 
and 7 becomes ill-conditioned. Even so, we can con- 
strain the parameter space by finding what combina- 
tions of parameters are consistent with the line pro- 
files. One question is if additional braking is necessary, 
or if a moderation of Vion by stretching known quanti- 
ties within their estimated uncertainties is sufficient to 
explain the profile. This should put a lower limit on 
the extra braking needed, and by extension put an in- 
dependent lower limit on the C II content in the disk. 
We get a lower reasonable limit for the Vion of Na I by 
assuming /3 330 (which is 1.5 cr below its estimated 
value) and an ionization rate that is 20 % higher than 
the nominal, resulting in vion = 2.5kms~^. In Fig. [4] 
the least square-sum of the residuals between the ob- 
served line profiles for Na I and Fe I, and their mod- 
els (where the system radial velocity again was a free 
parameter) are plotted as a function of 7/70, where 70 
corresponds to the /3 and F assumed above, an d a brak- 
ing ion gas density according to the model bv iZBWlOl 
The corresponding column density of C II is given on 
the upper horizontal axis of the figure. The smallest 
residuals are found for 7/70 = 0.5, corresponding to the 
C II column density A''cii = 1-3 x lO^^cm "^, which is 
2 X the model column density of IZBWlOl and 6 x the 
iVcii = 2.0to;4 X 10"cm-2 reported bv IROB061 The 
best-fit Na I velocity profile resulting from increased 
braking is shown in Fig. [2] Because the braking is pro- 
portional to the velocity, only the high-velocity tail of 
the profile is significantly affected by the C II density 
increase. Similarly, Fe I never gets accelerated to high 
enough velocities for braking to be important, so the Fe I 
velocity profile does not change with the increased C II 
density. 

The implicated overabundance of C is on the order of 
100 X compared to cosmic abundances of other detected 
elements. Perhaps this reflects the intr insic abundance 
of the gas production mechanism; e.g., iKarmann et al.l 
(j2003l ) find carbon rich FEBs essential to explain the 
high-velocity FEB spectral features. Another possibil- 
ity is that the current abundances are more strongly 
dependent on the gas removal mechanism. Even with 
an efficient ion braking mechanism, any specific parti- 
cle is neutral some fraction of the time. If the radia- 
tion pressure on the particles in their neutral state over- 
comes gravity, they will be accelerated outwards until 
they get ionized (and braked again). This will effec- 
tively remove elements at a rate proportional to their 
ionization speed Uion and neutral fraction. On the other 
hand, elements that do not experience sufficient radia- 
tion pressure in their neutral phase (i.e., H, He, C, N, O, 
F, Ne, and Ar, see Table 1 of FBW06) arc thus expected 
to pile up, thereby increasing their relative abundance 
compared to their production abundance. Unfortunately, 
the elements that experience little radiation pressure are 
also the elements that are difficult to observe, so only C 
and O of that group have so far been detected (,ROB061 , 
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although H has observed upper l i mits (jFreudling et alj 
Il995t iLecavelier des Etangs et all I2001D . Surprisingly, 
absorption line observation s by FU SE favor a cosmic 
abundance of disk O I gas (IROBOGT ). Again, this con- 
clusion is sensitively dependent on the assumed (unre- 
solved) line profile, as can be seen in Fig. 1 of theROJBOj 
online supplementary material: reducing the broadening 
parameter from the best-fit value of ~14kms~^ to a few 
kms^^ results in a much higher, but less constrained, 
column density. A broadening parameter as large as 
14kms^^ for the stable disk gas is unexpected in itself 
considering other observed absorption lines, and is likely 
caused by multiple unresolved velocity components from 
FEBs (as seen for other elements). 

One way to confirm the overabundance of C II around 
/3Pic is to observe its thermal 157.7 /xm emission, as 
the emission is predicted to be roughly proportional to 
the C gas mass (ZBWIO). The Infrared Space Obser- 
vatory (ISO) observed (3 Pic with the LWS instrument 
and found a 4 tr emissi on feature with an integrated flux 
of 10^13 gj.gs-icm-2 ()Kamp et al.ll2003[ ). This is 15x 
more emission than predicted from a gas disk with 20 x 
overabundant C (ZBWLQ), and would imply a C abun- 
dance of >100x over cosmic values (compared to other 
observed elements in the disk). Recently, the infrared 
space telescope Herschel with the PACS instrument ob- 
served both C II 157.7 /im and O I 63.2 /im emission, con- 
firming the ISO detection and finding that O I probably 
also is overabundant. The data are still being analyzed, 
however, as e.g. the spectrophotometric calibration is not 
settled yet. A publication presenting the Herschel obser- 
vations is in preparation by the key project Stellar Disk 
Evolution team lead by G. Olofsson. 



The processes that produce the C-rich gas disk around 
/3 Pic must reasonably be present also in other debris disk 
systems. The frequency and efficiency of the process is 
not known, but t he absence of detecta ble amounts of gas 
around AUMic (jRoberge et al.l l2005f). another edge-on 
debris disk system in the /3 Pic moving group, shows that 
the process is not ubiquitous. A model to quantitatively 
investigate the balance of gas production and removal 
mechanisms in debris disks will be the subject for future 
studies. 

6. CONCLUDING SUMMARY 



The prediction made by the braking model of'FBWd6|, 
that Fe I and Na I should show different radial veloci- 
ties, is confirmed. This shows that, to a first order ap- 
proximation, neutral elements are not braked by the gas, 
while the ions must be. Looking more closely at the line 
profiles, we see an absence of a Na I high- velocity popula- 
tion, expected from free acceleration. Our interpretation 
is that the high-velocity Na atoms must be braked. A 
factor of 2-5 X more C II gas than previously estimated 
is required for C II to be the braking agent; the implica- 
tion is that the gas disk around /3 Pic is even more C rich 
than previously thought, up to 100 x cosmic abundance 
relative to other detected elements. 
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